home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / studentbase.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  3KB  |  130 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. #define TWOVRPI 0.636619772367581343
  15. #define HALF_PI 1.5707963268
  16. #define EPSILON .000001
  17.  
  18. /* CACM Algorithm 395, by G. W. Hill */
  19.  
  20. void studentbase(x, df, cdf)
  21.         double *x, *df, *cdf;
  22. {
  23.   double t, y, b, a, z, j, n;
  24.  
  25.   n = *df;
  26.   z = 1.0;
  27.   t = (*x) * (*x);
  28.   y = t / n;
  29.   b = 1.0 + y;
  30.  
  31.   if (n > floor(n) || (n >= 20 && t < n) || (n > 20)) {
  32.     if (n < 2.0 && n != 1.0) {
  33.       /* beta integral aproximation for small df */
  34.       double  da = 0.5, db = 0.5 * n, dx, dp;
  35.       int ia = 0, ib = floor(db);
  36.   
  37.       dx = db / (db + da * t);
  38.       betabase(&dx, &db, &da, &ib, &ia, &dp);
  39.       *cdf = (*x >= 0) ? 1.0 - .5 * dp : .5 * dp;
  40.     }
  41.     else {
  42.       /* asymptotic series for large or non-integer df */
  43.       if(y > EPSILON) y = log(b);
  44.       a = n - 0.5;
  45.       b = 48.0 * a * a;
  46.       y = a * y;
  47.       y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5)
  48.             / (0.8 * y * y + 100.0 + b) + y + 3.0) / b + 1.0) * sqrt(y);
  49.  
  50.       y = -1.0 * y;
  51.       normbase(&y, cdf);
  52.       if (*x > 0.0) *cdf = 1.0 - *cdf;
  53.     }
  54.   }
  55.   else {
  56.     /* nested summation of cosine series */
  57.     if (n < 20 && t < 4.0) {
  58.       a = y = sqrt(y);
  59.       if(n == 1.0) a = 0.0;
  60.     }
  61.     else {
  62.       a = sqrt(b);
  63.       y = a * n;
  64.       for(j = 2; fabs(a - z) > EPSILON; j += 2.0) {
  65.         z = a;
  66.         y = (y * (j - 1)) / (b * j);
  67.         a = a + y / (n + j);
  68.       }
  69.       n += 2.0;
  70.       z = y = 0.0;
  71.       a = -a;
  72.     }
  73.     for(n = n - 2.0; n > 1.0; n -= 2.0) a = ((n - 1.0) / (b * n)) * a + y;
  74.     a = (fabs(n) < EPSILON) ? a/sqrt(b) : TWOVRPI * (atan(y) + a / b);
  75.     *cdf = z - a;
  76.     if(*x > 0.0) *cdf = 1.0 - 0.5 * *cdf;
  77.     else *cdf = 0.5 * *cdf;
  78.   }
  79. }
  80.  
  81. /* CACM Algorithm 396, by G. W. Hill */
  82.  
  83. double ppstudent(pp, n, ifault)
  84.      double pp, n;
  85.      int *ifault;     
  86. {
  87.   double sq, p, a, b, c, d, x, y;
  88.  
  89.   /* convert to double upper tailed probability */
  90.   p = (pp < 0.5) ? 2.0 * pp : 2.0 * (1.0 - pp);
  91.  
  92.   if (n <= 3.0) {
  93.     if (n == 1) sq = tan(HALF_PI * (1.0 - p));
  94.     else if (n == 2.0) sq = sqrt(2.0 / (p * (2.0 - p)) - 2.0);
  95.     else {
  96.       sq = ppbeta(p, 0.5 * n, 0.5, ifault);
  97.       if (sq != 0.0) sq = sqrt((n / sq) - n);
  98.     }
  99.   }
  100.   else {
  101.     a = 1.0 / (n - 0.5);
  102.     b = 48.0 / (a * a);
  103.     c = ((20700.0 * a / b - 98.0) * a - 16) * a + 96.36;
  104.     d = ((94.5 / (b + c) - 3.0) / b + 1.0) * sqrt(a * HALF_PI) * n;
  105.     x = d * p;
  106.     y = pow(x, 2.0 / n);
  107.     if (y > 0.05 + a) {
  108.       /* asymptotic inverse expansion about normal */
  109.       x = ppnd(0.5 * p, ifault);
  110.       y = x * x;
  111.       if (n < 5) c = c + 0.3 * (n - 4.5) * (x + 0.6);
  112.       c = (((0.05 * d * x - 5.0) * x - 7.0) * x - 2.0) * x + b + c;
  113.       y = (((((0.4 * y + 6.3) * y + 36.0) * y + 94.5) / c - y - 3.0) / b
  114.            + 1.0) * x;
  115.       y = a * y * y;
  116.       y = (y > .002) ? exp(y) - 1.0 : 0.5 * y * y + y;
  117.     }
  118.     else {
  119.       y = ((1.0 / (((n + 6.0) / (n * y) - 0.089 * d - 0.822)
  120.                    * (n + 2.0) * 3.0) + 0.5 / (n + 4.0)) * y - 1.0)
  121.         * (n + 1.0) / (n + 2.0) + 1.0 / y;
  122.     }
  123.     sq = sqrt(n * y);
  124.   }
  125.  
  126.   /* decode based on tail */
  127.   if (pp < 0.5) sq = -sq;
  128.   return(sq);
  129. }
  130.